#AUTHOR: Pierpaolo Da Fieno
#LICENSE: LGPL


from math import *

#Physical Constants

G=9.80665
PI=3.141593

#m^2/s
NI_15=0.118900E-05

#kg/m^3
RHO=1025

#kg/m^2
SIGMA_STD=5e6

#Controllable Pitch Propeller
CPP=False
CPP_COEFF=1.02


#
# rblades: Number of blades (Z)
# pitch2d: Pitch/Diameter ratio (P/D)
# expand: Expanded Area of propeller (Ae)
# advance: Advance coefficient (J)
# reynolds: Reynolds number (Re)



# Wageningen Coefficients

def kt(rblades,pitch2d,expand,advance):
    t1 = +0.00880496 \
    -0.204554*advance \
    +0.166351*pitch2d \
    +0.158114*pitch2d**2 \
    -0.147581*advance**2*expand \
    -0.481497*advance*pitch2d*expand \
    +0.415437*pitch2d**2*expand \
    +0.0144043*rblades \
    -0.0530054*advance**2*rblades \
    +0.0143481*pitch2d*rblades
    t2 = +0.0606826*advance*pitch2d*rblades \
    -0.0125894*expand*rblades \
    +0.0109689*advance*expand*rblades \
    -0.133698*pitch2d**3 \
    +0.00638407*pitch2d**6 \
    -0.00132718*advance**2*pitch2d**6 \
    +0.168496*advance**3*expand \
    -0.0507214*expand**2 \
    +0.0854559*advance**2*expand**2 \
    -0.0504475*advance**3*expand**2
    t3 = +0.010465*advance*pitch2d**6*expand**2 \
    -0.00648272*advance**2*pitch2d**6*expand**2 \
    -0.00841728*pitch2d**3*rblades \
    +0.0168424*advance*pitch2d**3*rblades \
    -0.00102296*advance**3*pitch2d**3*rblades \
    -0.0317791*pitch2d**3*expand*rblades \
    +0.018604*advance*expand**2*rblades \
    -0.00410798*pitch2d**2*expand**2*rblades \
    -0.000606848*rblades**2 \
    -0.0049819*advance*rblades**2
    t4 = +0.0025983*advance**2*rblades**2 \
    -0.000560528*advance**3*rblades**2 \
    -0.00163652*advance*pitch2d**2*rblades**2 \
    -0.000328787*advance*pitch2d**6*rblades**2 \
    +0.000116502*advance**2*pitch2d**6*rblades**2 \
    +0.000690904*expand*rblades**2 \
    +0.00421749*pitch2d**3*expand*rblades**2 \
    +0.0000565229*advance**3*pitch2d**6*expand*rblades**2 \
    -0.00146564*pitch2d**3*expand**2*rblades**2
    return t1+t2+t3+t4

def kq(rblades,pitch2d,expand,advance,cpp=CPP):
    t1 = +0.00379368 \
    +0.00886523*advance**2 \
    -0.032241*advance*pitch2d \
    +0.00344778*pitch2d**2 \
    -0.0408811*pitch2d*expand \
    -0.108009*advance*pitch2d*expand \
    -0.0885381*advance**2*pitch2d*expand \
    +0.188561*pitch2d**2*expand \
    -0.00370871*advance*rblades \
    +0.00513696*pitch2d*rblades
    t2 = +0.0209449*advance*pitch2d*rblades \
    +0.00474319*advance**2*pitch2d*rblades \
    -0.00723408*advance**2*expand*rblades \
    +0.00438388*advance*pitch2d*expand*rblades \
    -0.0269403*pitch2d**2*expand*rblades \
    +0.0558082*advance**3*expand \
    +0.0161886*pitch2d**3*expand \
    +0.00318086*advance*pitch2d**3*expand \
    +0.015896*expand**2 \
    +0.0471729*advance*expand**2
    t3 = +0.0196283*advance**3*expand**2 \
    -0.0502782*pitch2d*expand**2 \
    -0.030055*advance**3*pitch2d*expand**2 \
    +0.0417122*advance**2*pitch2d**2*expand**2 \
    -0.0397722*pitch2d**3*expand**2 \
    -0.00350024*pitch2d**6*expand**2 \
    -0.0106854*advance**3*rblades \
    +0.00110903*advance**3*pitch2d**3*rblades \
    -0.000313912*pitch2d**6*rblades \
    +0.0035985*advance**3*expand*rblades
    t4 = -0.00142121*pitch2d**6*expand*rblades \
    -0.00383637*advance*expand**2*rblades \
    +0.0126803*pitch2d**2*expand**2*rblades \
    -0.00318278*advance**2*pitch2d**3*expand**2*rblades \
    +0.00334268*pitch2d**6*expand**2*rblades \
    -0.00183491*advance*pitch2d*rblades**2 \
    +0.000112451*advance**3*pitch2d**2*rblades**2 \
    -0.0000297228*advance**3*pitch2d**6*rblades**2 \
    +0.000269551*advance*expand*rblades**2 \
    +0.00083265*advance**2*expand*rblades**2
    t5 = +0.00155334*pitch2d**2*expand*rblades**2 \
    +0.000302683*pitch2d**6*expand*rblades**2 \
    -0.0001843*expand**2*rblades**2 \
    -0.000425399*pitch2d**3*expand**2*rblades**2 \
    +0.0000869243*advance**3*pitch2d**3*expand**2*rblades**2 \
    -0.0004659*pitch2d**6*expand**2*rblades**2 \
    +0.0000554194*advance*pitch2d**6*expand**2*rblades**2
    if cpp:
        return CPP_COEFF*(t1+t2+t3+t4+t5)
    else:
        return t1+t2+t3+t4+t5
        
def delkt(rblades,pitch2d,expand,advance,reynolds):
    if reynolds<=2e6:
        delkt=0
    else:
        d1 = log10(reynolds)-0.301
        d2 = d1**2
        delkt = 0.000353485 \
        -0.00333758*expand*advance**2 \
        -0.00478125*expand*pitch2d*advance \
        +0.000257792*d2*expand*advance**2 \
        +0.0000643192*d1*pitch2d**6*advance**2 \
        -0.0000110636*d2*pitch2d**6*advance**2 \
        -0.0000276305*d2*rblades*expand*advance**2 \
        +0.0000954*d1*rblades*expand*pitch2d*advance \
        +0.32049E-5*d1*rblades**2*expand*pitch2d**3*advance
    return delkt

def delkq(rblades,pitch2d,expand,advance,reynolds):
    if reynolds<=2e6:
        delkq=0
    else:
        d1 = log10(reynolds)-0.301
        d2 = d1**2
        delkq = -0.000591412 \
        +0.00696898*pitch2d \
        -0.0000666654*rblades*pitch2d**6 \
        +0.0160818*expand**2 \
        -0.000938091*d1*pitch2d \
        -0.00059593*d1*pitch2d**2 \
        +0.0000782099*d2*pitch2d**2 \
        +0.52199E-5*d1*rblades*expand*advance**2 \
        -0.88528E-6*d2*rblades*expand*pitch2d*advance \
        +0.0000230171*d1*rblades*pitch2d**6 \
        -0.184341E-5*d2*rblades*pitch2d**6 \
        -0.00400252*d1*expand**2 \
        +0.000220915*d2*expand**2
    return delkq

def eta0(rblades,pitch2d,expand,advance):
    kq=kq(rblades,pitch2d,expand,advance)
    if kq == 0:
        return 0
    else:
        return advance/twopi*kt(rblades,pitch2d,expand,advance)/kq

def eta0rey(rblades,pitch2d,expand,advance,reynolds):
    kq=kq(rblades,pitch2d,expand,advance)+delkq(rblades,pitch2d,expand,advance,reynolds)
    if kq == 0:
        return 0
    else:
        return advance/twopi*(kt(rblades,pitch2d,expand,advance)+delkt(rblades,pitch2d,expand,advance,reynolds))/kq


#TODO:definitions

def t_sigma(rblades,rpm,diameter,pitch2d,powershaft,sigma=SIGMA_STD):
    powershaft_cv=powershaft/.7355
    sigma_kg_cm2=sigma/1E4
    result=( 40.5 *( 1- 0.474 *pitch2d ) * powershaft_cv / rblades / rpm / diameter**3 / ( sigma_kg_cm2 +0.059 * rpm**2 * diameter**2 ) )**( 1./3 )
    result*=.21
    return diameter * ( 0.0092 + result )

#c: chord of blade

def c(rblades,expand,diameter):
    return 2.073*expand/rblades*diameter

#t2c: Thickness to Chord Ratio

def t2c(rblades,expand):
    return (0.0185-0.00125*rblades)*rblades/2.073/expand

#t2c: Minimun Thickness to Chord Ratio for blade strength

def t2c_sigma(rblades,rpm,diameter,pitch2d,powershaft,expand,sigma=SIGMA_STD):
    result=t_sigma(rblades,rpm,diameter,pitch2d,powershaft,sigma)
    result/=c(rblades,expand, diameter)
    return result

#t2c_opt: Thickness to Chord Optimum ratio for cavitation

def t2c_opt(cavitation_no):
    return 0.3*cavitation_no-0.012

#Re: Thickness to Chord Optimum ratio for cavitation

def Re_eq(t2c,t2c_sigma,reynolds):
    return exp(4.6052+sqrt((1+2*t2c)/(1+2*t2c_sigma))*(log(reynolds)-4.6052))
    
def Re_eqJ(t2c,t2cs,c, speed_advance, advance,ni=NI_15):
    return Re_eq(t2c,t2cs,ReJ(c, speed_advance,advance,ni=ni))

def Re(c, diameter, rpm, speed_advance, ni=NI_15):
    return c/ni * sqrt( speed_advance**2 + (.75*PI*rpm/60*diameter)**2 )
    
def ReJ(c, speed_advance,advance,ni=NI_15):
    return c/ni * speed_advance/advance*sqrt( advance**2 + (.75*PI)**2 )

def ReJN(c, diameter, rpm, advance,ni=NI_15):
    return c/ni * sqrt( (rpm/60*advance*diameter)**2 + (.75*PI*rpm/60*diameter)**2 )

#TODO: definitions

def cavitation_no(diameter, rpm, speed_advance, water_height, density=RHO, r=0.7):
    vr2=speed_advance**2+(PI*rpm/60*diameter*r)**2
    return (99537+RHO*G*water_height)/(.5*RHO*vr2)
    
def tao_cavitation(cavitation_no):
    return 0.0095692+.4623*cavitation_no-0.19409*cavitation_no**2
    
def Tmax_cavitation(diameter, pitch2d, rpm, speed_advance, water_height, density=RHO, r=0.7):
    vr2=speed_advance**2+(PI*rpm/60*diameter*r)**2
    A=PI*diameter**2/4*(1.067-0.229*pitch2d)
    sigma_c=cavitation_no(diameter, rpm, speed_advance, water_height, density, r=r)
    return 2*tao_cavitation(sigma_c)*A*vr2

def tao(T,expand,diameter,pitch2d,rpm,speed_advance,density=RHO, r=0.7):
    vr2=speed_advance**2+(PI*rpm/60*diameter*r)**2
    A=PI*diameter**2/4*(1.067-0.229*pitch2d)
    return T/.5/density/vr2/A

